1 Set up

Install necessary packages:

install.packages("tidyverse", dependencies = TRUE) 
install.packages("raster", dependencies = TRUE) 
install.packages("sf", dependencies = TRUE) 
install.packages("landscapetools", dependencies = TRUE)
install.packages("landscapemetrics", dependencies = TRUE)

2 Prepare data

2.1 Quick visualisation

Import only the RGB color bands as individual RasterLayer objects:

library(raster)

#blue
b2 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B2.tif')

#green
b3 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B3.tif')

#red
b4 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B4.tif')

Combine the RasterLayer objects and visualise the satellite image:

landsatRGB <- stack(b4, b3, b2) #order is impt

plotRGB(landsatRGB, 
        stretch = "lin") #scale the values (try using "hist" also)

Landsat-8 true color composite (RGB). Source: U.S. Geological Survey.

2.2 Import data

Import all 5 bands from the satellite data as a RasterStack object named landsat:

filenames <- paste0('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B', 1:5, ".tif")

landsat <- stack(filenames)

#rename bands
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR')

Check coordinate reference system of landsat:

crs(landsat)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 48N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",105,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16048]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

2.3 Crop data

Import polygon of city boundaries as sgshp and check if their coordinate reference systems match:

library(sf)

sgshp <- st_read("data/master-plan-2014-region-boundary-web-shp/MP14_REGION_WEB_PL.shp")

Check coordinate reference system of sgshp:

crs(sgshp)
## [1] "PROJCRS[\"SVY21\",\n    BASEGEOGCRS[\"SVY21[WGS84]\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6326]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",1.36666666666667,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",103.833333333333,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",1,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",28001.642,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",38744.572,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"

Transform sgshp to the match the coordinate reference system of the landsat:

sgshp <- st_transform(sgshp, crs = crs(landsat))

Crop landsat according to city boundaries:

landsat <- crop(landsat, sgshp) #crop to rectangle

landsat <- mask(landsat, sgshp) #mask values according to shape of sgshp

Plot the cropped image using only the RGB bands:

landsatRGB <- subset(landsat, c(4,3,2)) #Red, Green, Blue

plotRGB(landsatRGB,
        stretch = "lin")

Landsat-8 true color composite (USGS, 2018) cropped to city boundaries (URA, 2014)

3 Classify land cover

3.1 Calculate NDVI

Create a function that calcuates the Normalized Difference Vegetation Index (NDVI) for each pixel:

ndvi <- function(x, y) {
    (x - y) / (x + y)
  }

Apply function to the NIR and Red bands of landsat

landsatNDVI <- overlay(landsat[[5]], landsat[[4]], 
                       fun = ndvi)

Limit the range of values to be from -1 to 1:

landsatNDVI <- reclassify(landsatNDVI, c(-Inf, -1, -1)) # <-1 becomes -1

landsatNDVI <- reclassify(landsatNDVI, c(1, Inf, 1)) # >1 becomes 1

3.2 Visualise NDVI

Map out the NDVI values:

plot(landsatNDVI, 
     col = rev(terrain.colors(10)), 
     main = "Landsat 8 NDVI",
     axes = FALSE, box = FALSE)

Plot histogram of NDVI values:

hist(landsatNDVI,
     main = "Distribution of NDVI values", xlab = "NDVI", 
     xlim = c(-1, 1), breaks = 100, yaxt = 'n')
abline(v=0.2, col="red", lwd=2)
abline(v=0, col="red", lwd=2)

3.3 Define NDVI threshold

Set 0.2 as the threshold; reclassify values below this threshold to NA:

landsatGreen <- reclassify(landsatNDVI, c(-1, 0.2, NA)) #-1 to 0.2 becomes NA

Plot values of NDVI larger than 0.2

plot(landsatGreen, 
     main = 'Vegetation cover',
     col = "darkgreen", 
     axes = FALSE, box = FALSE, legend = FALSE)

3.4 Classify using NDVI

Create a matrix to be used as an argument in the reclassify() function:

reclass_m <- matrix(c(-Inf, 0, 1, #water
                      0, 0.2, 2, #urban
                      0.2, Inf, 3), #veg
                    ncol = 3, byrow = TRUE)
reclass_m
##      [,1] [,2] [,3]
## [1,] -Inf  0.0    1
## [2,]  0.0  0.2    2
## [3,]  0.2  Inf    3

Classify land cover using the defined threshold values:

landsatCover <- reclassify(landsatNDVI, reclass_m)

Plot the land cover classes:

3.5 Save raster

Save the reclassified raster landsatcover in the GeoTiff format:

writeRaster(landsatCover, 
            filename = "clean_data/landsat_landcover.tif", 
            overwrite = TRUE)

4 Landscape metrics

4.1 Quick visualisation

library(landscapemetrics)
library(landscapetools)

#landsatCover <- raster('clean/landsat_landcover.tif') #reload saved raster

show_landscape(landsatCover, discrete = TRUE)

Check if the raster data is in the right format:

check_landscape(landsatCover)
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         3  ✔

4.2 Patch-level

E.g. Area of each patch in the landscape:

lsm_p_area(landsatCover)
## # A tibble: 12,551 × 6
##    layer level class    id metric value
##    <int> <chr> <int> <int> <chr>  <dbl>
##  1     1 patch     1     1 area   14.6 
##  2     1 patch     1     2 area    0.09
##  3     1 patch     1     3 area   27.3 
##  4     1 patch     1     4 area    0.09
##  5     1 patch     1     5 area    0.18
##  6     1 patch     1     6 area    0.36
##  7     1 patch     1     7 area    0.09
##  8     1 patch     1     8 area    0.09
##  9     1 patch     1     9 area    0.09
## 10     1 patch     1    10 area    0.09
## # ℹ 12,541 more rows

4.3 Class-level

E.g. For each class, the total area of all patches:

lsm_c_ca(landsatCover)
## # A tibble: 3 × 6
##   layer level class    id metric  value
##   <int> <chr> <int> <int> <chr>   <dbl>
## 1     1 class     1    NA ca      6737.
## 2     1 class     2    NA ca     28504.
## 3     1 class     3    NA ca     42930.

E.g. For each class, the average area of patches:

lsm_c_area_mn(landsatCover)
## # A tibble: 3 × 6
##   layer level class    id metric  value
##   <int> <chr> <int> <int> <chr>   <dbl>
## 1     1 class     1    NA area_mn  3.23
## 2     1 class     2    NA area_mn  4.52
## 3     1 class     3    NA area_mn 10.3

4.4 Landscape-level

E.g. Total area of the landscape (all three land cover classes):

lsm_l_ta(landsatCover)
## # A tibble: 1 × 6
##   layer level     class    id metric  value
##   <int> <chr>     <int> <int> <chr>   <dbl>
## 1     1 landscape    NA    NA ta     78170.

5 Credits

Spatial data used in this document:

Creative Commons Licence

© XP Song